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Abstract — MIMO processing techniques in fiber optical com- 
munications have been proposed as a promising approach to meet 
increasing demand for information throughput. In this context, 
the multiple channels correspond to the multiple modes and/or 
multiple cores in the fiber. The lack of back-scattering necessitates 
the modeling of the transmission coefficients between modes as 
elements of a unitary matrix. Also, due to the scattering between 
modes the channel is modeled as a random Haar unitary matrix 
between M transmitting and N receiving modes. In this paper, we 
apply a large-deviations approach from random matrix theory to 
obtain the outage capacity for Haar matrices in the low outage 
limit, which is appropriate for fiber-optical communications. 
This methodology is based on the Coulomb gas method for 
the eigenvalues of a matrix developed in statistical physics. By 
comparing our analytic results to simulations, we see that, despite 
the fact that this method is nominally valid for large number of 
modes, our method is quite accurate even for small to modest 
number of channels. 

Index Terms — Optical fiber transmission, MIMO, outage ca- 
pacity, random matrix theory 

I. Introduction 

THE ongoing exponential growth in wire-line data traffic 
is primarily driven by high-bandwidth digital applica- 
tions, such as video-on-demand, cloud computing and tele- 
presence. As a result, it is expected that the currently de- 
ployed infrastructure will soon reach its limits, leading to 
the so-called "capacity crunch" |[T|. To counter this trend, 
scientists have been working towards exhausting all available 
degrees of freedom of fiber-optical transmission, including 
the bandwidth (through WDM modulation), available power 
(subject to power constraints imposed by non-linearities), and 
polarization diversity |j2|. One additional possibility to increase 
throughput is spatial modulation, which would allow multiple 
transmission within the same fiber or fiber bundle. This can 
be achieved by designing multi-mode (MMF) and/or multi- 
core fibers (MCF). While such fibers may be constructed to 
have small cross talk between the propagating modes ||3J, in 
the typical situation cross-talk between fiber modes cannot be 
neglected, for reasons related to increased length of the fiber 
segment |j4J, as well as the twist and the bending of the fiber 
||5), ||6). In addition, the increase of the number of modes 
necessarily will bring the modes closer in space, resulting to 
higher cross-talk, which can lead to the power being spread 
equally between modes of the fiber Q. As a result, strong 
coupling between cores needs to be addressed with a different 
approach. 




Fig. 1 . Optical SDM using M parallel transmission paths in MCF. There is 
crosstalking between adjacent cores. 
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Fig. 2. Light beam scattering paradigm in MMF. 



Recently, Q, |[8), Q it has been proposed to use techniques 
developed in the context of wireless communications between 
multiple transmitting and receiving antennas (MIMO). As it 
is well known in wireless multi-antenna systems the crosstalk 
is significant and yet sophisticated techniques developed there 
can be used to reach large information throughputs. Of course, 
the application of MIMO methodology to optical fiber multi- 
core systems has certain differences which need to be ad- 
dressed. First, while the channel matrix for wireless MIMO is 
typically assumed to be complex Gaussian, this is not usually 
the case for fibers. It has been suggested instead that the 
channel matrix should be a subset of a unitary matrix ||2], 
Second, due to the existence of non-linearities at high 
powers, one should specifically have in mind low to moderate 
powers per channel. Hence, the qualitative behavior described 
in the context of large powers (usually called diversity multi- 
plexing tradeoff in the literature | |1 IJ ) is not valid. Finally, the 
practical metric for the performance is not the ergodic mutual 
information, but, rather, instead the outage capacity at very 
low outage (e.g. 10^^) |2|, due to the fact that feedback from 
the receiver to the transmitter is almost always impossible. 

In this paper, we analyze the outage capacity of an optical 
MIMO channel. As mentioned above, this is the relevant infor- 
mation transmission metric for fiber-optical coupled multi-core 
channels when the crosstalk varies over time slowly enough, 
so as not to be able to code around it. We obtain analytical 
expressions, which are valid technically in the limit of large 
channel numbers, but also works well over smaller channel 
numbers. It is particularly suited to obtain outage mutual 
information for very low outages with finite SNR. Essentially, 
it amounts to calculating the rate function of the logarithm 
of the average moment generating function of the mutual 
information. The methodology we use is based on the so-called 
Coulomb gas approach which was developed in the physics 
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literature in the context of random matrix theory p2) in the 
60's. It is quite intuitive because it interprets the eigenvalues 
as point charges on a line repelling each other logarithmically. 
The Coulomb gas method has seen recently a renewed interest 
in its use to obtain large deviations results for random matrix 
problems [13] , [14J and also in communications |15J , [16) . 
We will follow the basic steps discussed in more details in 

A. Outline 

In the next section we will define the system model, and 
show that that the appropriate channel matrix is a random Haar 
unitary matrix. In Section III we introduce the mathematical 
methodology of the Coulomb gas and provide our analytic 
results. In addition, we compare to numerical simulations. 
Finally, in Section IV we conclude. 

II. System Model 

In this paper we consider a single-segment A^-channel 
lossless optical fiber system, with Nt < N transmitting 
channels excited and Nr < N receiving channels coherently 
excited in the input (left) and output (right) side of the fiber 
The propagation through the fiber may be analyzed through 
its 2N X 2N scattering matrix given by 



4-T 



(1) 



Due to time-reversal symmetry this matrix is a complex unitary 
symmetric matrix (with S = S^) fl?) . The N x N blocks on 
the diagonal rt, represent the reflection matrices of the left 
ingoing and right ingoing channels. In the case of optical fiber 
transmission we are interested only in transmission from the 
left to right. Also proper manipulation of the input through 
the use of optical isolators results to no backscattering, i.e. 
rt = r,. ~ 0. Assuming strong scattering between right- 
moving channels, as discussed above, we may model t as 
a complex Haar-distributed matrix with t^^t = tt^ = I^r. 
We are particularly interested in a segment of this matrix 
corresponding to the Nt columns and the Nr rows, which 
are coupled to the transceiver. We denote this Nt x Nr matrix 
U and without loss of generality we assume it is the upper 
left corner of t. As a result, the corresponding MIMO channel 
for this system reads 



y = Ux 



(2) 



with coherent detection and channel state information only at 
the receiver | |T8| , p9) . x, y and z are the Nt x 1 input, the 
Nr X 1 output signal vectors and the iV^ x 1 unit variance 
noise vector, respectively, all assumed for simplicity to be 
complex Gaussian. We also assume no differential delays 
between channels, which effectively leads to frequency flat 
fading |j2). We also assume no mode-dependent loss. As a 
result, the mutual information can be expressed as 

In{V) = logdet(/ + pU^U) 

Nt 
fe=l 



(3) 



where "log" is the natural logarithm, p is the average total 
signal-to-noise ratio, are the eigenvalues of the matrix U^U 
and we assume for concreteness Nt < Nr. 



III. Methodology 

We are now able to define the main problem we address 
here, namely the calculation of 



Poutir) = Prob{lN < Ntr) 

= Eu[e{Ntr-lN{V))] 

We will first analyze the density of r i.e. 

Pir) = Poutir) 



(4) 
(5) 



(6) 



Due to space constraints, we will focus on the case Nt + 
Nr < N and expand on other parameter regions in a longer 
version. It is useful to define /3 = Nr/Nt > 1 and no = (N — 
Nt — Nr)/Nt The first step is to express the joint distribution 
of eigenvalues of U^U as derived initially in po) and more 
recently in this context pO| 

PA(Ai...AArJcx II |A„~A„p[]Af*-'^'-l(l-A,)^« 



iAre 



-N^E(\) 



(7) 



where An is a normalization constant and E{X) is an energy 
function of the eigenvalues A^. E{\) represents the potential 
energy of Nt unit charges bound to move on a the unit interval 
X e (0, 1), while repelling each other and from the boundaries 
logarithmically. It is reasonable to assume that when N is large 
then they will coalesce to a smooth density p{x), which will 
be such that the energy E{X) = £[p] will be minimum. Since 
the above equation is in fact a probability distribution, this 
will correspond to the most probable eigenvalue distribution. 
In this continuum limit we can rewrite the minimum energy 
as 



Sq = sup inf Cq [c, p] 

c P 

Cq = —uq / p{x) log(l — x)dx 



(8) 



— 1) J p{x)log{x)dx 
p{x)p{y) log \x - y\dydx 
p(x)dx — 1 

where we have added a Lagrange multiplier c to ensure that 
p{x) is properly normalized, while implicitly we assume that 
p{x) is continuous in x g (0, 1). The minimization of Co with 
respect to p can be done by taking its functional derivative 
and setting it to zero. More technical details about this can be 
found in |15) . Following this methodology, we can show that 
the local minimum of Cq is unique in the space of continuous 
i^+'^-integrable functions of p{x), where x G (0, 1). 
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As we will see from the analysis later on, the density 
function minimizing £o is 



Pq{x) 



y/{x - a){b - x) 
2ttx{1 — x) 



(9) 



where a, 6 = (VrTno±\/9(noT^)^/("-o + l+/3) obtained 
using other methods in pO| , pT) . We now move on to 
incorporate the constraint on the rate r by introducing a new 
Lagrange multiplier so that 



£{r) = supinf £[c, 

c,k P 



(10) 



C — Cq — k (^j p{x) log(l + px)dx ~ r 

The minimum of the above function is also unique since we 
can extract £q from £(r) maximizing over C while keeping 
k = |15|. The extra term above, proportional to k, has 
the intuitive explanation of an additional force acting on the 
charges-eigenvalues towards the left or the right, depending on 
the sign of k. Taking the functional derivative with respect to 
p{x) and setting to zero and differentiating with respect to x, 
provides us with the following integral equation representing 
a balance of forces on the charge located at x. 



2V 



p{x') 



dx' 



no 
1 - x 



kp 



1 



px 



(11) 



where V represents the Cauchy principal value of the integral 
and a,h € (0, 1) the limits of support of p. Following 
Tricomi's theorem |15|, [22 1 this integral equation may be 
solved to give a unique solution. For the values /? > 1 and 
no > we get 



p{x) 



yj{x~ a)(h-x) 
27r(l + px) 



1) 



a){l-b) 



(12) 



with the additional constraint (stemming from the continuity 
constraints p{a) = p{b) = 0) 



no 



/3-1 



kp 



v/(l-a)(l-6) y/{l + pa){l + pb) 



(13) 



The parameters a, b, k can be evaluated uniquely from the 
above equation in addition to the normalization constraint 

/ dxp{x) — 1 



no + /3 + l + fc = 



/3-1 
\/ab 



fc(l + p) 

v/(l + P«)(H 



pb) 



(14) 



and the rate constraint J dxp{x) log(l + px) 



log- 



-G 



-G 



A 



fcA 



a 



2^{z + a){z + b) 
f z 1 — a 



G 



a + z a + z 



A ' 
a + z 
A ' 



A 

1 - a 
A 



A 

(/3-l)A 
2^/ab 



G 



A 



(15) 



where A = 6 — a and z = -.The G function can be seen in 

p 

Appendix A. We may now integrate over p{x) and obtain an 
expression for £{t) as follows 



£ir) = - 



log A 



(no + /3 + 1) 
fcA 



2zy/{z + a){z + b 
13-1 



no 
2 



{Giw,y)~G{w,-<f>)) 



+G(0,2/)-G(0 
(/?-l)A 



2Vab 



no 
2 



(g (0, 4>)-G (0, -w 



{G{w,w)-G{w,<l>)) 



+G{0,w) - G (0,-0) 
fc 



2 

no 



- log 1 



2 log(l-a)-^^^loga 



(16) 



where cj) ^ w ^ -^,w = w + l,y = = 

y + I- £0 can be evaluated, e.g. from the above by setting 
fc = 0. 



A. Probability Distributions P{r) and Pout{f) 

We may now obtain P{r) from the above to leading 



exponential order, by properly normalizing e 



-Nt£{r) 



usmg 



the fact that due to exponential nature in N the bulk of the 
distribution will be centered around the minimum of £{r) with 
respect to r which occurs at fc = 0, which corresponds to 
Gaussian approximation. Hence 



P{r)^Nt- 



y^2nVe 



(17) 



where Verg = l/f"(7'erg) IS the variance at the peak of 
the distribution, and r^rg is the solution of (15 1 for fc = 
corresponding to the ergodic rate. It should be noted that r^rg 



1 23 1 and Verg | |24| have been calculated using other methods 
in the past. 

The outage probability may also be obtained in a similar 
fashion [15|. For r < re,.„ the outage probability is 



-N^l£i{r)-£o- 



Poutir) 



N\£'iir)\ 



(18) 



and for r > r^rg it is 



Poutir) « 1 



-N-[£,ir)-£o-l^]^f N\£[ir)\ 



(19) 



erg 



where £'i{r) and £f(r) are the first and second derivative of 
£i{r) with respect to r and Q{x) is the known erfc function. 
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N =10, N =6,p=10 
r t ^ 



f « 






LD Nj=2 




LD = 6 




Monte Carlo - 6 




Monte Carlo —2 




LD . -2 




Gaussian - 6 




Gaussian - 2 




Monte Carlo - 2 




Gaussian = -2 


1 
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(a) Outage probability curves for different values of A'o We observe that, generally,the Gaussian curves fail to follow the respective Monte Carlo, while 
the LD curves are closer to them. 



N =3,N, = 2,Ng = 2,p=10 



■LD 

■Gaussian 
■ Monte Carlo 




(b) For small values of p and Nr , Nt , Nq the discrepancy between the different methods is almost insignificant. 

N^ = 10, Nj = 6, p=200 



_LDN„ = 2 
.,.,LD Nu = 6 

Monte Carlo Njj = 2 

Monte Carlo N^^ = 6 

. . - Gaussian = 6 
Gaussian = 2 

TRT 




(c) At the extreme case of very high SNR, the discrepancy between the different methods is large. The TRT curve can successfully follow the LD and 
Monte Carlo curve only for low outage probabilities . 



Fig. 3. Simulation results 
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1) Numerical Simulations: We have performed a series of 
numerical simulations and have compared the Large Deviation 
(LD) approach to Gaussian approximation and Monte Carlo 
simulation. We start with the case of small p and for different 
values of Nt, and TVq. In the figures we plot the outage 
probability according to the LD, Gaussian and Monte Carlo 
approach, for r < r^rg- The curves are for different values 
of p and different realisations of MIMO systems. It is visible 
that the Pout from the Gaussian simulation induces significant 
deviation from the reality (Monte Carlo) at the "tails". As it 
is expected, for higher SNR, higher rates can be used, but 
an approach closer to realistic systems, "forces" us to select 
lower SNR (p). 

For large SNR (p = 200) we have inserted the 
Throughput Reliability Tradeoff (TRT) curve, . It is a piece- 
wise linear function of the outage probability logarithm 

log Pout (R) - c{k)R -g{k) log p (20) 
c(fc) = Nt+Nr-2k-l (21) 
g{k) = NtNr - k{k + 1) (22) 

where R = Ntr and is p is large and fclogp < R < (fc + 
l)logp 

IV. Conclusion 

In this paper we have used a large deviations approach first 
introduced in physics | ,12J to calculate the outage capacity for 
the optical MIMO channel in the limit of large channel num- 
bers. Our method is especially applicable for the tails of the 
distribution, which is relevant for low outage requirements due 
to the absence of feedback and finite SNR. Our analytical re- 
sults agree very well with numerical experiments. Additionally 
the method provides the distribution of eigenvalues constrained 
on the transmission rate and SNR. Although the channel 
assumptions taken here (random unitary without losses) are 
somewhat idealized, this result gives an analytic metric to 
compare with other more complicated channel models. 

Appendix 

The function G{x,y) for a; > and y > or y < — 1 is 
given by |[T5) 



(23) 



= -2sgn(?;)Vl2/(l + y)|log 
'y/\l + x\ + 



^\xil + y)\ + ^\y{l + x)\ 



+ |l + 2y|log 



V\^\ + V\y\ 
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